/**
 * mLibGaussianElimination.js
 * Created: Dec 2013
 * Source link: http://code.google.com/p/code-depository/source/browse/
 *
 * By Suresh Nalluri,
 * License: The MIT License (MIT)
 *          See LICENSE details
 *
 *
 *	Todo:

 * Requires:
 *	mLibIncludes.js

 * Functions used to solve an arbitrary sized linear system of equations
 * using Gaussian elimination.
 *
 * The equations are represented by a matrix which is stored in column major
 * order.
 * 
 * So we are solving Ax = b; Where A is an nxn matrix & b is an n-dim vector
 *
 *-------------------------------------------------------
 **/
if (typeof mLIB == 'undefined') function mLIB() {};

( function(module) {
    var abs = module.Includes.abs,
        isZero = module.Includes.isZero;
    
    module.gauss = function() {};  // declare the sub-module where these functions will be kept
    
    var friends = {
        hasSolution: function(A, b, n) {
            // Solving Ax = b, where A is an nxn matrix and b is nx1.
            //
            // return true or false depending on whether system is solvable or not
            // the actual solution is returned within the incoming vector b
            // note that the incoming matrix A representing the equations, gets 
            // clobbered.
            
            // basic algorithm:
            //
            // set the augmented matrix M is thus: [A | b]
            // using elementary row operations, we need to arrive at the following row echelon
            // composition, consisting of 1's in the diagonal
            //      +-                    -+
            //      | 1 v v ..... v | c1   |
            //      | 0 1 v v ....v | c2   |
            //      | ...           | .    |
            //      | ...           | .    |
            //      | 0 0 0 ... 0 1 | cn-1 |
            //      +-                    -+
            // 
            //
            // for each pivot row 0, 1, 2,.., n-1; the diagonal element is the pivot element.
            // determine the appropriate row from the matrix M, ie. the one who's diagonal
            // element has the largest absolute value, then swap it with the current pivot
            // row.  Once swapped, divide the now pivot row by its diagonal value, ie the now
            // pivot element; this will set the pivot element to 1
            //
            // for each row below the pivot row, we need to ensure the element e in the
            // column, which contains the above diagnol is set to 0, ie zero out the pivot column
            // elements below the pivot row: to do this, multiply the pivot row by e and subtract the row
            //
            // repeat above for next pivot row
            // note for the above steps the pivot row index is the same as the pivot column index
            //
            // now do backward substitution
            //
            //
            
            var iRow, pivot = 0, tmp, laVal, laRow, iCol, val, pivotRecip;
            
            // the elementary row ops first
            for (; pivot < n; pivot++) {
                /* get the row (laRow) with largest absolute element in the current column
                 */
                laRow = pivot; // the current row
                laVal = abs(A[laRow + n*pivot]);
                for (iRow = pivot+1; iRow < n; iRow++) {
                    val = abs(A[iRow + n*pivot]);
                    (val > laVal) && (laVal = val, laRow = iRow);
                }
                
                // if laVal is zero then we have a singular matrix .. no solutions exists
                if (isZero(laVal)) return false;
                
                /* swap the rows if laVal is not in the current row.
                 */
                if ( laRow != pivot ) {
                    for (iCol = 0; iCol < n; iCol++) {
                        tmp = A[pivot + n*iCol];
                        A[pivot + n*iCol] = A[laRow + n*iCol];
                        A[laRow + n*iCol] = tmp;
                    }
                    // don't forget to swap the rows in b too!!
                    tmp = b[pivot];
                    b[pivot] = b[laRow];
                    b[laRow] = tmp;
                }
                
                /* divide the current row by the pivot element
                 */
                pivotRecip = 1.0 / A[pivot + n*pivot];      // remember pivot element is always the diagonal
                for (iCol = 0; iCol < n; iCol++) {
                    A[pivot + n*iCol] *= pivotRecip;
                }
                b[pivot] *= pivotRecip;
                
                // let's set pivot to 1 anyway so we have a clean 1
                A[pivot + n*pivot] = 1.0;
                
                /* zero out all the elements in pivot column, below the current (pivot) row.
                 */
                for (iRow = pivot+1; iRow < n; iRow++) {
                    val = A[iRow + n*pivot];
                    
                    // multiply the pivot row with val and subtract from this row
                    // we start with iCol = pivot; since all elements of the row before this column
                    // were zeroed from previous iterations.
                    for (iCol = pivot; iCol < n; iCol++) {
                        A[iRow + n*iCol] -= val*A[pivot + n*iCol];
                    }
                    b[iRow] -= val*b[pivot];
                }
            }
            
            // backward substitution;  we'll use a do while to illustrate backwards substitution
            var iRow = n-1;  // start at the bottom row; note we have b[n-1] already by default
            do {
                iRow--;
                for (iCol = iRow+1; iCol < n; iCol++) {
                    b[iRow] -= A[iRow +n*iCol]*b[iCol];
                }
            } while ( iRow > 0 );
            
            return true;
            
        },  // hasSolution()
        
        
        invertMatrix: function(A, n) {
            // invert given matrix A (given in col-major order).  We put the inversion within A
            // as we reduce the matrix to the reduced row echelon form
            // At the end of the gauss elim steps, we need to unswap (any row swaps) but we swap the
            // corresponding columns
            var iRow, pivot = 0, tmp, laVal, laRow, iCol, val, pivotRecip,
                swap = new Array[3];

            // the elementary row ops first
            for (; pivot < n; pivot++) {
                /* get the row (laRow) with largest absolute element in the current column
                 */
                laRow = pivot; // the current row
                laVal = abs(A[laRow + n*pivot]);
                for (iRow = pivot+1; iRow < n; iRow++) {
                    val = abs(A[iRow + n*pivot]);
                    (val > laVal) && (laVal = val, laRow = iRow);
                }
                
                // if laVal is zero then we have a singular matrix .. no solutions exists
                if (isZero(laVal)) return false;
                
                /* swap the rows if laVal is not in the current row.
                 */
                swap[pivot] = laRow;   // note laRow could be same as pivot
                if ( laRow != pivot ) {
                    for (iCol = 0; iCol < n; iCol++) {
                        tmp = A[pivot + n*iCol];
                        A[pivot + n*iCol] = A[laRow + n*iCol];
                        A[laRow + n*iCol] = tmp;
                    }
                }
                
                /* divide the current row by the pivot element
                 */
                // pivotRecip
                pivotRecip = 1.0 / A[pivot + n*pivot];      // remember pivot element is always the diagonal
                for (iCol = 0; iCol < n; iCol++) {
                    A[pivot + n*iCol] *= val;
                }
                
                // set pivot to pivotRecip.. this is where we build the inverse in-line
                // we wont be changing this bit again as next iteration moves us to the next row
                A[pivot + n*pivot] = pivotRecip;
                
                /* zero out all the other elements in pivot column, below and above the current (pivot) row.
                 */
                for (iRow = 0; iRow < n; iRow++) {
                    if ( iRow == pivot ) continue;
                    
                    val = A[iRow + n*pivot];
                    
                    // clear the pivot column (building the inverse in-line)
                    // it will eventually get set to -val*pivotRecip
                    A[iRow + n*pivot] = 0.0;
                    
                    // multiply the pivot row with val and subtract from this row
                    // this should set the 
                    for (iCol = 0; iCol < n; iCol++) {
                        A[iRow + n*iCol] -= val*A[pivot + n*iCol];
                    }
                }
            }
            
            // now we undo the swaps if any, in reverse order
            iCol = n;  // start at the bottom row; note we have b[n-1] already by default
            
            do {
                iCol--;
                if (swap[iCol] != iCol) {
                    // note we swap the corresponding columns of the matrix
                    for (iRow = 0; iRow < n; iRow++) {
                        tmp = A[ iRow + n*iCol];
                        A[ iRow + n*iCol] = A[iRow + n*swap[iCol]];
                        A[iRow + n*swap[iCol]] = tmp;
                    }
                }
            } while ( iCol > 0 );
            
            return true;
            
        },  // invertMatrix
        
        determinant = function(M, n) {
            // some properties of determinants that will make for a better understanding of the method here
            //  1. det(AB) = det(A)det(B)
            //  2. if U is upper triangular matrix then det(U) = product of diagonal elements.  and if all
            //     diagonals are 1 (as is the case with a row echelon form matrix, this gives det(U) = 1
            //  3. using 1 and 2: SUppose that S1, S2, ..Sn are matrices that represent elementary row ops
            //     performed on A to reduce it to an upper triangular matrix.
            //          ie. [Sn][Sn-1]...[S1]A = U then we have the following
            //              det([Sn][Sn-1]..[S1]A) = det(U) = 1
            //              => det([Sn][Sn-1]..[S1]).det(A) = 1
            //              => det(A) = 1 / ( det(Sn)det(Sn-1)...det(Sn) )
            //
            //  now as we transform A to U, the following holds for the intermediate determinants:
            //  i)   multiplying a row by scalar k:         multiplies determinant by k
            //  ii)  Adding multiples of row to another row: no effect
            //  iii) Swapping rows:                         changes the sign of the determinant
            //
            var A = M.slice(0), // make a copy, we mangle this rather than destroy A
                det = 1.0,      // set initial value of determinant
                pivot = 0,
                iRow, iCol,
                pivotRecip, tmp, val;

            // lets reduce the matrix to the upper traingular form
            for (; pivot < n; pivot++) {
                /* get the row (laRow) with largest absolute element in the current column
                 */
                laRow = pivot; // the current row
                laVal = abs(A[laRow + n*pivot]);
                for (iRow = pivot+1; iRow < n; iRow++) {
                    val = abs(A[iRow + n*pivot]);
                    (val > laVal) && (laVal = val, laRow = iRow);
                }
                
                // if laVal is zero then we have a singular matrix .. determinant = 0
                if (isZero(laVal)) return 0.0;
                
                /* swap the rows if laVal is not in the current row.
                 */
                if ( laRow != pivot ) {
                    for (iCol = 0; iCol < n; iCol++) {
                        tmp = A[pivot + n*iCol];
                        A[pivot + n*iCol] = A[laRow + n*iCol];
                        A[laRow + n*iCol] = tmp;
                    }
                    det = -det; // row swap so negate determinant
               }
                
                /* divide the current row by the pivot element
                 */
                pivotRecip = 1.0 / A[pivot + n*pivot];      // remember pivot element is always the diagonal
                for (iCol = 0; iCol < n; iCol++) {
                    A[pivot + n*iCol] *= pivotRecip;
                }
                det *= pivotRecip;
                
                // set pivot to 1 anyway so we have a clean 1
                A[pivot + n*pivot] = 1.0;
                
                /* zero out all the elements in pivot column, below the current (pivot) row.
                 */
                for (iRow = pivot+1; iRow < n; iRow++) {
                    val = A[iRow + n*pivot];
                    
                    // multiply the pivot row with val and subtract from this row
                    // we start with iCol = pivot; since all elements of the row before this column
                    // were zeroed from previous iterations.
                    for (iCol = pivot; iCol < n; iCol++) {
                        A[iRow + n*iCol] -= val*A[pivot + n*iCol];
                    }
                }
            }

            return 1.0 / det;
            
        } // determinant()
        
    };  // friends{}
    
    
    for (var f in friends) {
        friends.hasOwnProperty(f) && (module.gauss[f] = friends[f]);
    };

    
} )(mLIB);